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Abstract 

We provide a pedagogical introduction to numerical linked-cluster expansions (NLCEs). We sketch the algorithm for 
generic Hamiltonians that only connect nearest-neighbor sites in a finite cluster with open boundary conditions. We 
then compare results for a specific model, the Heisenberg model, in each order of the NLCE with the ones for the 
finite cluster calculated directly by means of full exact diagonalization. We discuss how to reduce the computational 
cost of the NLCE calculations by taking into account symmetries and topologies of the linked clusters. Finally, we 
generalize the algorithm to the thermodynamic limit, and discuss several numerical resummation techniques that can 
be used to accelerate the convergence of the series. 

Keywords: Linked-cluster expansions, Exact diagonalization, Spin systems, Lattice models 



1. Introduction 

1.1. High-temperature expansions 

Calculating exact finite-temperature properties of 
quantum lattice models can be very challenging. One 
general approach to achieve that goal is to devise se- 
ries expansions for the lattice model in question in 
the thermodynamic limit |Q], 0] ■ A common example 
of these series expansions are high-temperature expan- 
sions (HTEs), in which the partition function X and 
other extensive properties of the system are expanded in 
powers of the inverse temperature f3 = (A^T) -1 (in what 
follows we set the Boltzman constant ks to unity). For 
example, consider the Ising Hamiltonian with nearest- 
neighbor interactions: 



parameter in the expansion: 
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where / is the strength of the interaction, (..) denotes 
nearest neighbors, and <x, (= +1) is the Ising spin on 
site i. The partition function can be written as 



(2) 



where the sum runs over all possible configurations of 
the spins. We define K = f3J, which serves as a small 
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If we associate (c^cry)' to an /-fold bond between sites 
; and j, a typical term in the above expansion can be 
represented graphically as the lattice with various num- 
ber of lines (including no line) connecting every two 
nearest-neighbor sites. Therefore, one can write the par- 
tition function in terms of contributions from all pos- 
sible graphs, up to a certain size, that can be embed- 
ded on the lattice. (In the Ising model, the fact that 
o~ — + 1 makes the calculations much easier than in 
quantum models such as the Heisenberg model.) The 
order in K that each graph contributes to is equal to the 
number of bonds it has (see Ref. J2l for more details 
on this derivation). Equation (O implies that the series 
converges only at high temperatures, above or of the or- 
der of /. In what follows, we will see how this type of 
expansion is related to linked-cluster expansions. 

1.2. Low-temperature expansions 

Similar to HTEs, low-temperature expansions (LTEs) 
can be devised to describe properties of a system with 
discrete excited states close to its ground state, i.e., for 
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P — > oo. In this approach, the partition function is writ- 
ten as 



Z 



-m 
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where Eq (£„) is the ground state (« th excited state). If 
there is an energy increment, e, satisfying E„ - Eq - 
m„e with m„ being integer for all n, any thermodynamic 
property can be expressed as an expansion in powers 
of e~^. Then, cluster expansions similar to the ones 
discussed above for HTEs follow J2l- 

1.3. Linked-cluster expansions 

The idea behind linked-cluster expansions (LCEs) Q 
l^t] is that for any extensive property P of a lattice model 
(such as the logarithm of the partition function or the 
internal energy) one can compute its value per lattice 
site P(£)/N in the thermodynamic limit in terms of a 
sum over contributions from all clusters c that can be 
embedded on the lattice: 



P(£)/N = Y j L(c)xW p (c), 



(5) 



where L(c) is the multiplicity of c, namely, the number 
of ways per site in which cluster c can be embedded on 
the lattice, and Wp(c) is the weight of that cluster for the 
property P. Wp(c) is defined according to the inclusion- 
exclusion principle: 



W, 



P (c) = P(c) - Yj Wits), 



where 



P(c) 



Tr[P(c)e-^] 
Tr [e-M] 



(6) 
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is the property calculated for the finite cluster c and the 
sum on s runs over all sub-clusters of c. In Eq. (0, 
H c is the Hamiltonian of cluster c. One can check that 
the weight of a disconnected cluster vanishes because 
P(c) can be written as the sum of its parts (see Sec. 12.21 ). 
hence, the name linked-cluster expansions. 

The convergence of the series in Eq. ©, when all the 
terms are considered, is guaranteed by the definition of 
weights in Eq. ((6]). In fact, by swapping the place of 
P(c) and Wp(c) in Eq. (O, one can write the property of 
a finite cluster c as the sum of its weight and the weights 
of its sub-clusters. Taking the thermodynamic limit c — > 
X. brings one back to Eq. (01. However, in that limit, 
only a finite number of terms can be accounted for, and 
the series has to be truncated. 

Because of the inclusion-exclusion principle 
[Eq. ((§)], the weight of every cluster contains only 



the contribution to the property that results from 
correlations that involve all the sites in the cluster, and 
in a unique fashion in accord with its specific geometry. 
At low temperature, when correlations grow beyond the 
size of the largest clusters considered in the series, the 
results show a divergent behavior (due to the missing 
contributions of clusters in higher orders of the expan- 
sion). In most of the two-dimensional quantum models 
of interest, this occurs near or at zero temperature, 
e.g., for the nearest-neighbor antiferromagnetic (AF) 
Heisenberg model on a bipartite lattice. 

The clusters in the sum (0) are usually grouped to- 
gether based on common characteristics to represent 
different orders J3,|5t]. For instance, in the site expan- 
sion, where sites are used as building blocks to generate 
the clusters, the order of the expansion is determined by 
the number of sites of the largest clusters. All the clus- 
ters with n sites are said to belong to the n order. In 
LCEs, one has the freedom to devise an expansion (with 
a certain building block for generating the clusters in 
different orders) that best suites the particular model of 
interest. Some of these include the site, bond, or square 
expansions on the square lattice, and site, bond, or tri- 
angle expansions on the triangular and Kagome lattices, 
and so on @]. 

Despite the simple form of the LCE equations above, 
its computational implementation can be challenging, as 
one has to (i) generate all the linked clusters that can be 
embedded on the lattice, (ii) identify their symmetries 
and topologies to compute multiplicities [this step dra- 
matically reduces the computational effort as, for any 
given model, many different clusters have identical val- 
ues of P(c)], (iii) identify the sub-clusters of each clus- 
ter to calculate weights, and (iv) calculate the property 
of each cluster and perform the sums. LCEs are compu- 
tationally very demanding as the number of embedded 
clusters, and their sub-clusters, grow exponentially with 
increasing the order of the expansion. Below, we ex- 
plain all those steps in the context of an example (the 
site expansion on a finite square lattice). We should 
stress that the HTE explained above for the Ising model 
can be seen as a LCE in which the property for each 
cluster is calculated using a perturbative expansion of 
Eq. (0 in terms of K. 

1.4. Numerical linked-cluster expansions 

In this work, we present a pedagogical overview of 
the numerical linked-cluster expansions (NLCEs) intro- 
duced in Ref. 10] ■ NLCEs use the basis of LCEs, but 
employ exact diagonalization (ED), instead of pertur- 
bation theory as done in HTEs, to calculate the prop- 
erties of finite clusters in the series. This means that 
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the properties of each cluster are calculated to all or- 
ders in /?. The main advantage of NLCEs over HTEs 
is that, for models with short-range correlations, it is 
possible to access arbitrarily low temperatures because 
of the lack of a small parameter in the series. Further- 
more, for models in which correlations grow slowly as 
the temperature is lowered, NLCEs can converge well 
below the temperature at which HTEs diverge. In NL- 
CEs, as opposed to HTEs, the convergence temperature 
is controlled by the correlations in the model, and by the 
highest order in the series that can be calculated. 

The basics of NLCEs, and results for various spin and 
itinerant models in the square, triangular and kagome 
lattices, were presented in Refs. Ji], [5], [(J . Recent ap- 
plications of this method exploring properties of frus- 
trated magnetic systems can be found in Refs. 
NLCE studies of the Hubbard model in the square and 
honeycomb lattices were reported in Refs. llOlll llll2ll . 
How to deal with Hamiltonians and observables that 
break some of the symmetries of the underlying lat- 
tice was discussed in Refs. lf]~3L [Till . Finally, how to 
generalize NLCEs to calculate ground-state as well as 
low-temperature properties of lattice Hamiltonians with 
dimerized ground states was the subject of Ref. II 1511 . 
while direct ground state calculations in the thermody- 
namic limit were done in Ref. 01611 and dynamics were 
explored in Ref. fl7ll . 

Here, we introduce NLCEs for finite clusters, in or- 
der to show how they converge to the exact result with 
increasing the order in the expansion, and, later, discuss 
NLCEs in the thermodynamic limit. The exposition is 
organized as follows. In Sec. [2] we present the algorith- 
mic details and the numerical implementation of NLCE 
in two dimensions for a finite 4x4 lattice. In Sec. 12.71 
we report results of the expansion for the AF Heisenberg 
model on this cluster and compare them, in each order, 
to exact results that can be obtained by means of full ex- 
act diagonalization of the 4x4 lattice. In Sec. [3] we dis- 
cuss how to extend NLCEs to the thermodynamic limit, 
and compare NLCE results for the Heisenberg model 
to those from large-scale quantum Monte Carlo (QMC) 
simulations. In Sec. [4] we describe two resummation 
techniques that have been found to be widely applicable 
to accelerate the convergence of the NLCEs for thermo- 
dynamic properties of lattice models of interest. 

2. Implementation of NLCEs for Finite Systems 

In this section, we sketch the NLCE algorithm for 
an arbitrary Hamiltonian that only connects nearest- 
neighbor sites on a N = 16 site square lattice with open 
boundary conditions, which is shown in Fig.Q] We have 



Figure 1 : The 4x4 finite lattice with open boundary conditions used 
as an example in the derivation of the NLCE. 



chosen this small lattice because it allows us to carry out 
the NLCE to all orders in the site expansion. It also al- 
lows us to compare the properties in each order of the 
expansion to exact results calculated directly by ED of 
the entire 16-site system. 

2.1. Embedded clusters 

There exists ( N ) — clusters in the « th order on 

\n) n\(N—n)\ 

a lattice with N sites. To identify them in a computer, 
one can number the lattice sites in an arbitrary fashion. 
Then, construct an array of size N whose ; th element is a 
binary number representing site number i on the lattice. 
A as an element of this array indicates that the corre- 
sponding site is not part of the generated cluster whereas 
a 1 indicates that it is part of the cluster. Therefore, all 



Table 1: Total number of embedded clusters (second column), number 
of linked-clusters (third column), number of linked-clusters that are 
not related by point-group symmetries (fourth column), and number 
of linked-clusters that are topologically distinct (fifth column), in each 
order of the site expansion for a 1 6-site lattice with nearest-neighbor 
interactions (Fig.fTJ. 
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clusters in the « th order can be generated by exploring 
all possible configurations of 1 and as elements of the 
above array, while keeping the total number of nonzero 
elements fixed at n. 

In the second column of Table [1] we list for our 
example of the 16-site lattice in Fig. Q] Once we have 
all the site clusters, we need to connect the sites present 
on them with bonds. This is done following the Hamil- 
tonian of interest. For models with nearest-neighbor in- 
teractions, the case of interest here, bonds are added to 
every pair of nearest-neighbor sites fl. Hence, each site 
in our clusters is connected by bonds with up to 4 other 
sites. Examples of clusters within the site expansion are 
given in Fig. [2] 

2.2. Connected clusters 

Many of the site clusters generated in the step above 
contain disconnected parts. Those clusters should be 
discarded because their weight Wp(c), as given by 
Eq. ©, is zero. This can be easily seen if one assumes 
that cluster c consists of two disconnected sub-clusters 
c\ and c%. Then, since P(c) is extensive, it can be written 
as the sum of the properties of the two sub-clusters: 



P(c) = P( Cl ) + P(c 2 ). 



(8) 



But, in addition to all of their subclusters, c\ and co 
themselves are among the subclusters of c. Therefore, 



starting from one of the sites in the cluster, we add bonds 
between that site and all of its neighboring sites, pro- 
vided that they are part of the cluster. The same process 
is repeated for each of the new sites and this continues 
until there are no more options for adding bonds. The 
resulting cluster is then compared to the original one. 
If they are not identical, the cluster must have discon- 
nected parts and is discarded. 

In the third column in Table Q] we show the number 
of clusters in each order of the site expansion after the 
disconnected ones have been removed. Figure|2l depicts 
some of the linked clusters generated for our 16-site lat- 
tice. 

2.3. Point-group symmetries 

If the Hamiltonian preserves the symmetries of 
the underlying lattice, clusters related by point-group 
symmetries have the same weight for all observables. 
Therefore, they can be grouped together in order to 
avoid repeating the calculation of the weights for all of 
them. Depending on the lattice, one may have different 
number of point-group symmetries. For the square 
lattice, which is the case in our example, one has the 
following eight point group symmetries: 



W, 



p{c) = P^-^Wpis) 



Identity 

Rotation by 90° 
Rotation by 180° 
Rotation by 270° 



Reflection about x = 
Reflection about y — 
Reflection about x = y 
Reflection about x = —y 



= P(c)- W P ( Cl ) + J] Wp(s) 



W P (c 2 ) + J] Wp(s) 



P(c) - P( Cl ) - P(c 2 ) = 0. 



(9) 



One then needs to find a way to distinguish between 
connected and disconnected clusters. A cluster is con- 
nected if starting from any site one can reach any other 
site moving along the bonds (only nearest-neighbor 
bonds in our case). We have implemented this idea by 
reconstructing the cluster from scratch. In our codes, 



i L-rx n 



Figure 2: A sample of linked clusters that can be embedded on our 
finite 16-site lattice. Clusters 3 and 5 are related by point-group sym- 
metry and are topologically the same as cluster 4. 



Hence, at this point in the implementation of the NL- 
CEs, the goal is to identify and keep only those clus- 
ters that are not related by point-group symmetries. We 
achieve this through the following steps, which need to 
be followed independently for each order of the expan- 
sion: 

(i) Create a list of symmetrically distinct clusters and 
their multiplicity. This list begins with the first 
linked cluster in the order under consideration with 
multiplicity one. 

(ii) Take the next connected cluster and generate all 
clusters that are symmetrically related to it by ap- 
plying the symmetry operations mentioned above. 

(iii) For each cluster generated in (ii), check whether it 
is present in the list created in (i). This is achieved 
by finding out whether there is a translation vector 
that converts one cluster into the other. If yes, in- 
crease the multiplicity of the stored cluster that is a 
match by one, and go to (ii). 
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(iv) If none of the clusters generated in (ii) is in the list 
created in (i), store the connected cluster in the list, 
set its multiplicity to one, and go to (ii). 

This step significantly reduces the number of clusters 
that one has to consider. Among the clusters shown in 
Fig.|2j clusters 3 and 5 are related by a 90° rotation and 
only one of them needs to be stored. The fourth col- 
umn in Table [TJ shows the number of clusters that are 
not related by point-group symmetries in each order of 
the site expansion. 

2.4. Topological clusters 

Furthermore, even if some clusters are not related 
by point-group symmetries, their Hamiltonian may be 
identical. For example, in models with only nearest- 
neighbor terms, clusters 3 and 4 in Fig.|2]have the same 
Hamiltonian. We then say that these two clusters are 
topologically equivalent, and only one of them needs to 
be diagonalized. It is important to note that topologi- 
cally equivalent clusters share the same thermodynamic 
properties, but they may lead to different results for cor- 
relation functions. For example, the distance between 
the two extreme sites in clusters 3 and 4 in Fig. [2] is 
2a and V2a (a is the lattice spacing), respectively, and 
that difference needs to be taken into account when cal- 
culating correlation functions. Still, since the full exact 
diagonalization of each cluster is the most time consum- 
ing part of the NLCE calculations, the fact that one only 
needs to diagonalize topologically different clusters (or 
simply, topological clusters) reduces the computation 
time dramatically. 

At this point, we need to generate a list of all topolog- 
ical clusters. It can be easily verified that they share the 
same topologically distinct sub-clusters, i.e., only diag- 
onalizing topological clusters allows one to carry out the 
entire NLCE calculations. 

The identification of the cluster topologies can be 
done through adjacency matrices (M). An adjacency 
matrix contains all the information about the spatial re- 
lations between sites in the cluster. In the simplest form, 
an adjacency matrix of an n-site cluster for a model with 
only nearest-neighbors interaction can be a n x n ma- 
trix whose rows/columns represent different sites and 
the matrix elements, Mjj, are either 1, if sites i and j 
are nearest neighbors (are connected by a bond), or 
otherwise. Such matrix will be symmetric with zeros as 
diagonal elements. A generalized version of such ma- 
trices can be used for models with bond anisotropy or 
correlations beyond nearest neighbors. 

If two clusters of the same size are topologically 
equivalent, then there exists a permutation of sites that 



transforms the adjacency matrix of one onto the other. 
Therefore, similar to the approach taken in the previous 
step, we make a list of topologically distinct clusters for 
each order as follows: 

(i) Create a list of topological clusters and their mul- 
tiplicity. This list begins with the first symmetri- 
cally distinct cluster in the order under considera- 
tion with its multiplicity. 

(ii) Take the next symmetrically distinct cluster and 
construct its n \ possible adjacency matrices (given 
the n\ site permutations). 

(iii) For each adjacency matrix generated in (ii), check 
whether it coincides with the adjacency matrix of 
any of the clusters in the list created in (i). If yes, 
increase the multiplicity of the stored cluster that is 
a match by the multiplicity of the new cluster, and 
go to (ii). 

(iv) If none of the adjacency matrices generated in (ii) 
matches the one of a cluster in the list created in 
(i), store the newly found topological cluster in the 
list, set its multiplicity to the one it had as a sym- 
metrically distinct cluster, and go to (ii). 

In practice, this part of the algorithm can be very time 
consuming not only because the number of permuta- 
tions can be very large, but also because of the matrix 
comparisons. So, alternatively, one can generate a key 
(an individualized number) for each adjacency matrix 
and use that key for comparison purposes. To do that, 
one has to pick only one, out of the n\ adjacency ma- 
trices to represent a cluster. This can be accomplished 
by labeling sites according to a unique criterion for the 
adjacency matrix. For example, if we think of each ad- 
jacency matrix as a large binary number by attaching 
its columns/rows one after another, one can pick the 
permutation of site numbers that yields the largest (or 
smallest) binary number. In Fig. [3] we show one of the 
adjacency matrices of a four-site cluster and the same 
matrix after exchanging site labels 1 and 2. The latter, 
which yields the smallest binary number that consists of 
columns of the matrix, is the one that will represent this 
cluster and that will be stored. 

The recipe for generating the keys is arbitrary and, for 
large cluster sizes, one may not be able to find a recipe 
that guarantees a unique key for each cluster, i.e., mul- 
tiple adjacency matrices may end up sharing the same 
key. The following is an example of one such recipe for 
creating the keys: 

key -^.V/,,(/x. /) •'. (10) 

U 



5 



3 t 



M : 
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10 11 
10 
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key = 1113 



Figure 3: An example of the adjacency matrix that represents a 4-site cluster (left). Rows and columns are representatives of the sites on the cluster. 
A 1 (0) as the (i, j) element of the matrix indicates that sites (' and j are connected (not connected) by a nearest-neighbor bond. The matrix on the 
right is sorted by exchanging sites 1 and 2. A key is associated to the sorted matrix according to an arbitrary formula [here, we use Eq. 4 1 OH . 



If the key is not unique, then when the key of a new 
cluster matches the one of a stored topological cluster, 
one has to directly compare the adjacency matrices in 
order to determine whether the two clusters are topo- 
logically the same. We emphasize that the use of keys 
accelerates the process of comparing adjacency matri- 
ces tremendously and has been implemented in all our 
NLCEs. 

Up to this point, we have identified and stored the 
topological clusters of each size, their multiplicities, and 
their adjacency matrices with their respective keys. The 
number of topological clusters in each order of our ex- 
ample is shown in the last column in Table. Q] This ta- 
ble now makes apparent the significant reduction of the 
number of clusters that need to be fully diagonalized 
from the initial r^j. In the following step, we explain 
how to enumerate the sub-clusters of each topological 
cluster to be able to ultimately calculate its weight. 

We should stress that the approach described here to 
identify topologically distinct clusters works so long as 
the cluster sizes are small enough that all permutations 
of the introduced adjacency matrices can be explored. 
For larger cluster sizes, one needs to generate more so- 
phisticated adjacency matrices whose description is be- 
yond the scope of this introduction to NLCEs. 

2.5. Sub-clusters 

Given the description so far for finding topologi- 
cal clusters, identifying all subclusters of a topologi- 
cal cluster is a relatively straightforward task. Note that 
finding the topological clusters of our finite 16-site lat- 
tice can, in itself, be interpreted as finding all subclus- 
ters of a 4x4 cluster in the list of topological clusters of a 
larger lattice. Therefore, following the same procedure 
as in Sec. l2.1l and Sec. 12.21 and replacing the 4x4 system 
with any of the topological clusters, we first generate 
all possible connected subclusters. This time, there is 
no need for checking the point-group symmetries of the 
subclusters (Sec. 12.3b . This is because all clusters that 
are related by symmetry share the same topology and 
have the same adjacency matrix. So, after generating 



a new connected subcluster, it is sufficient to construct 
its adjacency matrix and compare it to those of clusters 
with the same size that are stored in the list of topolog- 
ical clusters. Since there will always be a match, one 
only needs to keep track of how many subclusters of 
a certain topology a topological cluster has. The latter 
can be achieved by considering a multiplicity matrix, Y, 
whose (i, j) element gives the number of times the y'th 
topological cluster can be embedded in the ;th topolog- 
ical cluster. 

The steps described from Sec. l2.ll to Sec. l2.5l need to 
be carried out just once for all Hamiltonians involving 
only nearest-neighbor terms in the square lattice. The 
table with all topological clusters and subclusters can 
be stored and used for different nearest neighbor Hamil- 
tonians and, for any given Hamiltonian, for different mi- 
croscopic parameters. 

2.6. Weights 

Now, we have all the necessary tools to account 
for the subcluster subtractions in order to calculate the 
weights in Eq. ([6j. The steps that follow need to be car- 
ried out every time that a new Hamiltonian or Hamilto- 
nian parameter is explored. One starts with the smallest 
cluster in the first order. That cluster has no subclusters 
and, therefore, the weight of the property P in the cluster 
is simply equal to the property P(l). Then, to compute 
the weight of the next topological cluster (a single bond 
in the second order of our site expansion, see Fig. 0, 
we have to subtract the weight of the cluster in the first 
order from its property, according to the multiplicities 
given by the matrix Y: 

Wp(l) = P(l) 

W P (2) = P(2)-Y 21 W P (l) = P(2)-Y 2l P(l) 

: (ID 

In the above equations, indices 1 and 2 refer to the clus- 
ter number, c, which is not necessarily the order number. 
As mentioned in Sec. 11.31 in NLCEs, the property P(c) 
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is calculated using ED. For instance, the internal energy 
is: 



E(c) = (H c ) = 



Zfj l£ ,,(c)ex P [-/fe„(c)] 
ex P[-/?e«(c)] 



(12) 



where H c is the Hamiltonian for the finite cluster c with 
eigenvalues s„(c), and M c is the size of the Hilbert space 
on that cluster. 

We define partial sums, S , to group together contri- 
butions from topological clusters with a common char- 
acteristic in one order of the NLCE. For example, in 
the site expansion, the most natural characteristic is the 
number of sites. Therefore, S„, which is the sum of the 
contributions to property P from all n-site topological 
clusters in the series (c„) reads 

S„ =Y J L{c n )W P {c n ). (13) 
The property in the mth order of NLCE is then 

m 

P m (£)/N = Y J S n . (14) 
11=1 

2. 7. Results for the AF Heisenberg model 

Now that we have sketched the site expansion for a 
16-site lattice and for generic Hamiltonians that con- 
nect only nearest-neighbor sites, let us consider a spe- 
cific example of one such Hamiltonians, namely, the AF 
Heisenberg model 



(15) 



where S, is the spin operator at site i. For simplicity, we 
set J = 1 so that, in what follows, all energies are given 
in units of /. 

In Fig. 3] we show results for the energy per site (a), 
the entropy per site (b), 



lnZ + 



and the specific heat per site (c), 

1 (H 2 ) - (H) 2 



C v 



N 



T 2 



(16) 



(17) 



on the 16-site lattice. The exact results (labeled "Ex- 
act" in Fig. |4]i were obtained by means of full diag- 
onalization of the Hamiltonian, and are compared to 
those obtained in different orders of the site expansion 
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Figure 4: (a) Energy, (b) entropy, and (c) specific heat per site vs 
temperature for the AF Heisenberg model on the 16-site cluster in 
Fig- n Exact results are obtained by full exact diagonalization of the 
16-site cluster with open boundary conditions and are compared to 
those in different orders of the site expansion explained throughout 
this work. 



[Eq. (fT4l il. Two things are apparent in those plots: (i) 
with increasing order, the NLCE results converge to the 
exact ones at lower temperatures, (ii) At any given or- 
der, quantities that can be represented by lower order 
derivatives of the partition function converge to lower 
temperature. Note the contrast between the results for 
the energy and C v . Interestingly, even in the 15th order, 
when only the contribution from the 16-site cluster is 
missing, the energy has a visible discrepancy with the 
exact result at temperatures as high as T = 0.3. 



3. Thermodynamic Limit 

The generalization of the NLCE presented above to 
the thermodynamic limit is straightforward. Only cer- 
tain steps in the algorithm change when dealing with 
the infinite-size system. In this section, we discuss 
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those changes and show NLCE results for the Heisen- 
berg model [Eq. (T5[ ] in the thermodynamic limit. 

First, for the generation of the clusters, we start with 
the smallest one (a site in the site expansion) and sys- 
tematically add more building blocks (sites in the site 
expansion). Hence, to generate clusters with n + 1 sites, 
we consider all symmetrically distinct clusters that have 
n sites and add a nearest-neighbor site to every site at 
the edge of the cluster. Note that this way of building 
site clusters works because of translational invariance in 
the infinite lattice, and it is not appropriate for the finite 
system in Fig. Q] where not all sites are equivalent. This 
approach not only guarantees the generation of all possi- 
ble clusters that can be embedded on the infinite lattice, 
but also produces only connected clusters, i.e., the step 
presented in Sec . 12.21 (examining the connectivity of the 
clusters) can be skipped. The second column in Table|2] 
show the number of connected clusters per site, with up 
to 17 sites, in the site expansion of the square lattice in 
the thermodynamic limit. 

For finite systems without translational symmetry, 
all the different embeddings of a particular cluster are 
automatically generated in the approach discussed in 
Sec. 12. II Thus, the multiplicity of a symmetrically dis- 
tinct cluster is calculated by counting the times it ap- 
pears in the process of generating all clusters. For the 
infinite system, because of translational symmetry, the 
multiplicity is simply equal to the number of point- 
group symmetries that transform a cluster to another 
cluster that is not related to the first one by a translation. 
Hence, a process similar to the one carried out for iden- 
tifying the symmetries of the clusters on a finite lattice 
needs to be implemented for the clusters of the infinite 
system. However, in the latter, the multiplicity of sym- 
metrically distinct clusters are determined immediately 
after building them by just examining their point-group 
symmetries. Moreover, if by adding sites to a cluster in 
the previous order one finds a new cluster that is related 
by a symmetry operation to one of the clusters stored in 
the new list, we simply dismiss it. The third column in 
Table|2]show the number of symmetrically distinct clus- 
ters per site, with up to 17 sites, in the site expansion of 
the square lattice in the thermodynamic limit. 

After finding all symmetrically distinct clusters, we 
follow exactly the same procedure described in Sees. 12.41 
and 12.51 to determine all topological clusters and their 
subclusters. For the site expansion in the square lattice, 
the last column in Table [2] show the number of topolog- 
ical clusters per site in the thermodynamic limit. We 
should stress that all the steps described above need to 
be carried out just once for all lattice Hamiltonians that 
only connect nearest-neighbor sites in the square lattice. 




0.6 
<*> 0.4 
0.2 




0.6 











■ yf 








O ;'• 





0.4 



0.2 



(O ' \\\ 1 




ooo \\ \ 




i '• 












'/ / \ 




■ a % 








, , , iii , i 





0.3 



Figure 5: (a) Energy, (b) entropy, and (c) specific heat per site of the 
AF Heisenberg model on the square lattice in the thermodynamic limit 
as a function of temperature. NLCE bare sums up to 15th order are 
compared with QMC results for a 256 X 256 lattice. 



Once the topological clusters and subclusters are de- 
termined, one can study any nearest-neighbor model of 
interest. In Fig. [5] we show results for the energy (a), the 
entropy (b), and the specific heat (c) per site for the AF 
Heisenberg model on the square lattice, and up to the 
15th order of the site expansion. NLCE results are com- 
pared to those obtained by means of QMC simulations, 
previously reported in Ref. lfl3ll . for a very large peri- 
odic system with 256 x 256 sites. Those plots exhibit 
qualitatively similar features to the ones in Fig. [4] for 
a finite cluster. Namely, as the order increases, NLCE 
results converge to lower temperature. In addition, the 
energy and the entropy can be seen to converge to lower 
temperature than the specific heat. We note that, only 
in the 15th order, there are 42,192 topological clusters 
that need to be diagonalized (see Table [2}. Because of 
the very large number of clusters, one can use an embar- 
rassingly parallel code that distributes groups of clusters 
to different processors so that one diagonalizes many of 



them at once in every order of the NLCEs. 

4. Resummation Algorithms 

As seen from the results in Fig.|5j the bare sums of the 
NLCE convergence down to a temperature that depends 
on the order up to which the expansion can be carried 
out and, on a more fundamental level, it depends on the 
build up of correlations in the system as the tempera- 
ture is lowered. The longer the correlations the larger 
the cluster sizes that need to be included in the series 
to achieve convergence. Fortunately, even if one cannot 
calculate more orders of the NLCE, because of the ex- 
ponential increase of the number of clusters and of the 
Hilbert space of each cluster, several numerical resum- 
mation algorithms have been developed that accelerate 
the convergence of NLCEs ^ . 

The goal of those algorithms is to estimate P{£.)/N = 
limm^ooPm in Eq. (fT4l > from a finite set \P m }. Resumma- 
tion techniques then provide results for the observables 
of interest in regions where the bare sums [Eq. ( TPfl il 
do not converge. Here we will focus our discussion on 
two such methods: Wynn's algorithm ll 811 and Euler's 
transformation lfl9ll . They have been shown to be ex- 
tremely useful in accelerating the convergence of NL- 
CEs for thermodynamic properties in several models of 




Table 2: Total number of linked-clusters (second column), number of 
linked-clusters that are not related by point-group symmetries (third 
column), and number of linked-clusters that are topologically distinct 
(fourth column) per site, up to the 17th order of the site expansion of 
the square lattice (in the thermodynamic limit) with nearest-neighbor 
interactions. 
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Figure 6: (a) Energy, (b) entropy, and (c) specific heat per site of the 
AF Heisenberg model on the square lattice as a function of tempera- 
ture. Here we show results of last two orders (thick lines are used for 
the last order and thin lines for the one to last order) for the NLCE 
bare sums (up to the 15th order), Wynn's algorithm with 6 cycles of 
improvement, and Euler's transformation (for the last 1 1 orders). Re- 
sults from a large-scale QMC and exact diagonalization of a 1 6-site 
cluster with periodic boundary conditions are also shown. 
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In Wynn's algorithm, one defines 

s n — 0, s^ n — P n , 

Jk) _ Jk-2) 1 

" " +1 Ae* _1) ' 



(18) 



where the differentiating operator A is applied to sub- 
scripts and is expressed as 



(19) 



(21) 

It is expected that the even entries e„ converge to 
P(-QIN, while the odd ones ej, 2 ' +1) usually diverge, 
where / is defined as the number of cycles of improve- 
ment. As an example, assume that the highest order in 
the NLCE that can be considered is m, i.e., the bare sum 
yields P,„ for the property in the last order. Now, assume 
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that we use the Wynn algorithm with I — 1 (one cycle of 
improvement) instead. In that case, 



1 - 



(20) 



where the second term in the right hand side is the first 
order correction to the result of the bare sum. Yet, it can 
often significantly improve the convergence. Note that, 
for every cycle of improvement, there will be two terms 
less in the new series determined by means of Eq. ([LSl l 
{P m — > s m 2l ). 

We have found that this algorithm is one of the best 
for accelerating the convergence of NLCEs. In Fig.|6ja), 
we show the same energy as in Fig.[3a) but include the 
results obtained after six cycles of improvement using 
Wynn's algorithm. One can see that the last two orders 
of the Wynn sum (red solid lines) agree with each other 
down to much lower temperatures (T ~ 0.2) in compar- 
ison to the bare sums. In addition, the last order (thick 
solid line) is in perfect agreement with the QMC results 
in the entire temperature region shown. Very similar 
results can be seen for the entropy in Fig. HJb). Fur- 
thermore, as depicted in Fig. |6jc), and unlike the bare 
sums, the Wynn algorithm can precisely capture the po- 
sition and the height of the peak in the specific heat of 
this model, and even converges to temperatures lower 
than those accessible to the QMC calculations depicted 
in the figure. 

Note that the maximum cluster sizes considered in 
NLCE shown there (15 sites) are far smaller than the 
size of the lattice in the QMC simulations (256 x 265). 
Hence, we see that a thorough exploration of all topolo- 
gies in NLCEs, up to those cluster sizes, completely 
eliminates finite-size effects. In addition, in combina- 
tion with resummation algorithms, one can compute 
observables up to quite low temperatures. NLCE re- 
sults can also be contrasted against those obtained from 
ED of a 16-site cluster with periodic boundary condi- 
tions. As seen in Fig. [6] the finite cluster results depart 
from QMC at temperatures higher than the convergence 
temperature of NLCE even with the bare sum. Using 
slightly larger clusters, which can still be solved by 
means of full exact diagonalization, does not improve 
this trend either |l3 ] . 

We should point out that the AF Heisenberg model on 
the square lattice is one of the most challenging models 
that one can attempt to solve with NLCEs. This is be- 
cause the antiferromagnetic correlations grow exponen- 
tially with decreasing temperature, and quickly exceed 
the linear size of the largest clusters that can be treated 
within ED. For this reason, frustrated magnetic systems, 



for which QMC techniques run into the infamous sign 
problem, and where correlation lengths are small and 
grow slowly with decreasing temperature, are ideal to 
be explored by means of NLCEs. 

Another very useful resummation technique for alter- 
nating series, i.e., series in which S„ alternates in sign, 
is the Euler transformation. Within this approach, S„ is 
replaced by u„ = (-l)"S„ and P(-Q/N is approximated 
by the following sum: 



«o — U\ + «2 + 



- + X yr A ' M <» (21) 

1=0 1 



where A is defined as forward differencing operator 
A u n — u n+ \ — u n , 



A u„ - u„+2 - 2m„ + i + u n , 

A 3 u n = m„ + 3 - 3u„+2 + 3m„ + i - u n , 



(22) 



and « - 1 is the number of terms for which a bare sum 
is performed before the Euler transformation sets in for 
higher order terms. The first few terms of the approxi- 
mation can be written as 



(23) 



It is evident from the results in Fig. [6] after imple- 
menting the Euler transformation for the last 1 1 terms 
(n = 5), that the latter algorithm also dramatically im- 
proves the convergence from that of the bare sum. 
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